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Micro-local,  non-linear,  resolution  analysis 
of  generalised  Radon  transform  inversions 
in  anisotropic  media 

Maarten  V.  de  Hoop  and  Norman  Bleistein 

Center  for  Wave  Phenomena,  Colorado  School  of  Mines,  Golden  CO  80401-1887,  USA. 

Abstract  The  resolution  analysis  of  generalized  Radon  Transform  (GRT)  inversion  of 
seismic  da>a  is  carried  out  in  general  anisotropic  media.  The  GRT  inversion  formula  is 
derived  form  the  ray-Bom  tqrproximation  of  the  wave  field  for  volume  scatterers.  However,  by 
considering  scattering  surfaces  in  the  resolution  analysis,  rather  than  parameter  perturbations, 
we  show  that  the  inversion  provides  a  reflectivity  map  and  reflection/transmission  coefficients 
as  functions  of  scattering  angles  and  azimuths.  Those  coefficients  can  be  subjected  to  any  type 
of  Amplitude  Versus  Angle  (AVA)  or  Anq)litude  Versus  Offset  (AVO)  analysis.  By  applying 
the  inversion  to  Kirchhoff  approximate  data  rather  than  Bom  approximate  data,  we  show  that 
the  ou^t  is  actually  linear  in  the  reflection  coefficients  and,  hence,  a  nonlinear  function  of  the 
change  in  medium  parametos  across  discontinuity  surfaces — the  reflectors  of  the  medium. 


1.  Introduction 

Current  acquisition  techniques  take  reflection-seismic  measurements  at  increasingly  larger 
scattering  angles  both  for  image  enhancement  and  improved  estimation  of  elastic  parameters. 
To  achieve  a  better  resolution,  however,  one  must  account  for  the  interplay  between 
heterogeneity  and  anisotropy  at  different  length  scales.  Also,  to  accommodate  inversion  of 
wide  angle  data,  one  has  to  take  account  of  the  nonlinear  dependence  of  the  reflected  amplitude 
on  changes  in  medium  parameters,  at  least  near  the  relevant  reflectors  in  the  configuration.  We 
use  the  generalised  Radon  transform  or  GRT  formulation  to  achieve  these  goals. 

Our  analysis  begins  with  the  ray-Bom  approximation  of  the  direct  scattering  problem, 
which  represents  the  scattered  field  as  a  volume  integral.  This  representation  is  recast  as  a  sum 
of  surface  integrals  over  the  discontinuity  surfaces  of  the  medium  parameters — ^the  reflectors 
in  the  subsiuface.  Then  the  inversion  of  the  linearised  scattering  problem  is  carried  out  The 
resolution  analysis  of  the  linearised  GRT  inversion  serves  as  the  basis  to  arrive  at  a  micro-local 
nonlinear  formulation;  its  range  of  validity  is  extended  by  applying  a  stationary  phase  analysis 
with  respect  to  migration  dip,  defined  as  the  normalised  gradient  of  total  travel  time  along  the 
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scattering  characteristic.  In  this  process,  carrying  out  the  GRT  inversion  on  Kirchhoff-type 
data,  an  explicit  adjustment  of  the  inversion  formula  is  found. 

By  identifying  the  data,  for  each  image  point,  in  common  scattering-angle/azimuth  gath¬ 
ers,  the  GRT  inversion  formula  can  be  modified  to  obtain  refiection/transmission  coefficients 
at  specular  ray  geometries.  Any  amplitude  versus  scattering  angle — ^AVA — analysis  can  then 
be  applied  to  those  coefficients  to  derive  information  about  the  medium  perturbation.  Several 
parametrizations  of  this  perturbation  have  been  developed  to  reveal  to  which  quantities  the 
coefficients  are  sensitive. 

The  idea  of  estimating  angle-dependent  reflectivity  has  been  developed  by  various  authors; 
see  for  example  Geoltrain  and  Chovet  [1]  and  Lumley  [2].  A  more  rigorous  discussion  of 
such  an  approach  can  be  found  in  Bleistein  et  al.  [3].  In  our  paper,  we  derive  a  GRT-based 
inversion  procedure  that  accomplishes  the  goal  of  constructing  full  reflection — and  possibly 
transmission — coefficients  as  functions  of  scattering  angle  and  azimuth,  in  closed  form. 

GRT-type  inversion  formulas  for  the  linearised  inverse  problem  have  been  developed  by 
Norton  and  Linzer  [4],  by  Beylkin  [5, 6, 7,  8],  by  Miller  et  al.  [9, 10],  by  Beylkin  et  al.  [11], 
and  by  Rakesh  [12]  for  the  acoustic  case.  The  extension  to  the  elastic  case  was  discussed  by 
Beylkin  and  Burridge  [13],  and  anisotropy  was  considered  by  De  Hoop  etal.  [14]  and  Spencer 
and  E>e  Hoop  [15].  Inversion  formulas  aiming  to  estimate  reflectivity  rather  than  the  medium 
perturbation  were  developed  by  Cohen  and  Bleistein  [16],  by  Bleistein  and  Cohen  [17],  and  by 
Bleistein  [18, 19].  This  paper  brings  the  two  approaches  together.  Discussion  of  the  numerical 
implementation  of  GRT  inversion  procedures  can  be  found  in  De  Hoop  and  Spencer  [20]. 

The  GRT  approach  is  a  highfrequency  approach  to  inversion.  There  are  two  equally  valid 
points  of  view  about  the  utility  of  this  method.  First,  for  full  bandwidth  data,  we  obtain  by  this 
method  only  the  most  singular  part  of  the  solution  to  the  inverse  problem,  since  it  is  this  part 
of  the  solution  that  is  tied  to  the  limit  of  frequency  approaching  infinity.  On  the  other  hand,  as  a 
practical  matter,  the  typical  bandwidth  of  the  seismic  inverse  problem  is  such  that  the  data  can 
be  viewed  as  highfrequency  data  for  most  of  the  length  and  time  scales  of  the  geophysical 
model.  Thus,  numerical  implementation  of  the  derived  inversion  formulas  produce  useful 
results  for  seismic  exploration. 

Throughout  the  paper  we  have  excluded  the  possibility  of  multipathing  of  the  rays 
connecting  the  image  point  to  a  source  and  a  receiver  in  the  predefined  background  medium. 
Certain  changes  in  the  formula’s  are  necessary  in  that  case,  which  we  will  investigate  in  a 
separate  paper. 
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2.  The  basic  equations 

2.1.  Notation 

First,  we  introduce  some  basic  notation.  Choose  coordinates  in  the  configuration  according  to 
X  =  (xi,  X2,  X3)  =  Cartesian  position  vector, 

8  =  {si,S2,  S3)  =  source  point, 

*•  =  (^i5  r2,  ra)  =  receiver  point, 
t  =  time. 

The  medium  is  described  by 

p(x)  =  density, 

=  elastic  stiffness  tensor, 
while  the  wavefield  is  described  by 

u(x,t)  =  (ui(x,t),U2(x,t),u3(x,t))  =  displacement  vector, 
and  generated  by  a  source  distribution  given  by 

f(x,t)  =  (fi(x,t),  f2(x,t),  f3(x,t))  =  body-force  source  density. 

In  the  remainder  of  the  paper,  we  will  employ  the  summation  convention. 

/ 

2.2.  Governing  wave  equation 

The  displacement  in  a  heterogeneous  anisotropic  medium  satisfies  the  elastodynamic  wave 
equation 

pd'^ui  -  dj{cijktdiuk)  =  /,• ,  (1) 

with  summation  over  repeated  lower  case  indices,  here  and  below.  Let 

G{x,x',t)  =  {Gij,{x,x\t))  (2) 

be  the  causal  Green’s  tensor,  which  satisfies  (cf.  Eq.(l)) 

pdtGip  -  dj(cijkedtGkp)  =  SipS{x  -  x')S{t) ,  Gip  =  0  for  t  <  0  .  (3) 

2.3.  Asymptotic  ray  theory 

Here,  we  summarize  the  formulation  of  anisotropic  ray  theory  for  the  evaluation  of  the  Green’s 
tensor  (see,  e.g.,  Kendall  etal  [21]).  Let 

N 


terms  smoother  in  t . 


(4) 
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In  this  equation,  the  arrival  time  and  the  associated  polarization  vector  satisfy 


(P  Sik  -  ef)  =  0  (at  all  *) ,  (5) 

which  implies  the  eikonal  equation 

det(p  Sik  -  cijktidtr)  (OjT))  =  0  (at  all  x) .  (6) 

The  polarization  vectors  are  assumed  to  be  normalised  so  that  =  1.  Define  the 

slowness  vector  7^^^  by 

=  (7) 

Then,  Eq.(6)  constrains  7  to  lie  on  the  sextic  surface  A{x)  given  by 

det(/?^,jt  -  Cijktjtlj)  =  0  •  (8) 


A{x)  Consists  of  three  sheets  >t^^)(x),  H  =  1,2,3,  each  of  which  is  a  closed  surface 
surrounding  the  origin.  An  individual  sheet  is  also  described  by  (cf.  Eq.(5)) 


2%  =  p-  ^iCijkeinj^k  =  0  .  (9) 

The  scalar  amplitudes  must  satisfy  the  transport  equation 

Sifewer’d"’  S/tW) = 0 ,  (10) 


where  N,  again,  indicates  the  mode  of  propagation,  that  is,  the  sheet  of  the  slowness  surface 
on  which  the  corresponding  slowness  vector  lies. 

TTie  characteristic  or  group  velocities  are  normal  to  >l(^)(x)  at  7^^)  and  satisfy 


=  1 


V  = 


7 

n^Q 


see  Eq.(9).  The  normal  or  phase  speeds  are  given  by 


v(^)  = 


1 


(11) 


(12) 


The  unit  phase  direction  follows  as 

^(JV)  ^  y{N)^(n) 


From  Eq.(ll)  it  follows  that 

yW  =  Ir*"'!  cosx*"* , 

where  is  the  angle  between  and  7^^!. 

The  amplitudes  can  be  expressed  in  terms  of  certain  Jacobians, 


|t,(x')|F(*) 


At:\p{x)p{x')MYI'^ 


with  M  = 


dx  dx 
dqi  ^  dq2 


A  — 
dqt  dq2 


(13) 


(14) 
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in  which  A,  M,  v  and  V  carry  the  superscript  (N).  Here,  (91,92)  parameterize  the  rays 
originating  from  the  source.  One  can  verify  that  the  dimension  of  A  is  [time]^  x  [mass]“S 
which  upon  multiplication  by  force,  with  dimensions  [mass]  x  [length]  x  [time]“^,  gives  the 
dimension  of  displacement,  i.e.,  [length]. 

2.4.  Source  and  receiver  Green ’s  functions 

In  the  integral  representation  for  the  scattered  field,  we  need  the  Green’s  functions  originating 
both  at  the  source  and  the  receiver  points.  Further,  the  gradient  of  total  travel  times  from  the 
source  to  a  scattering  point  to  the  receiver  are  required  in  preparation  of  the  GRT  inversion. 
We  introduce  these  functions  here. 

Set 


G{x,t)  =  G(x,8,t) ,  G{x,t)  =  G{r,x,i) . 

Employing  asymptotic  ray  theory  in  both  Green’s  functions,  we  introduce  the  notation 

s)  ,  A^^\x)  =  A^^\r,  x) 

in  the  case  of  scattering  from  incident  mode  N  to  outgoing  mode  M. 

According  to  Eq.(7),  the  slowness  vectors  at  a;  are  given  by 

a) ,  7<^)(a;)  =  Va;T<^)(r,  aj)  ; 
the  associated  phase  directions  are  given  by 


^(A/)  =  J. _ 


and  the  phase  speeds  (cf.  Eq.(12))  are  given  by 

1  1 


y(N)  _ 


V(M)  ^ 


(15) 


(16) 


(17) 


(18) 


(19) 


■  i¥"'l 

We  also  define  the  two-way  travel  time  and  its  gradient, 

T(^^)(r,  y,  s)  =  r(^)(y,  s)  +  r(^)(r,  y)  , 

r(^^)(r,  *,  a)  =  V*r(^^)(r,  aj,  a) 

From  Eq.(17)  we  see  that 

r<^^)(r,a;,a)  =  ^^^(x)  -|- ^^^(aj) . 

The  direction  of 

*'=  |r(iVM)|  ’ 

will  be  the  nugration  dip,  which  we  referred  to  in  the  introduction.  The  ray  geometry  and 
slowness  vectors  are  illustrated  in  figure  1;  the  associated  polarization  vectors  are  shown  in 
figure  2. 


(20) 


(21) 
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3.  Medium  description:  micro-local  perturbation 

Let  c{x)  denote  a  smoothly  varying  background  medium.  To  analyse  a  typical  geological 
setting,  we  consider  the  following,  micro-local,  representation, 

c<*)(*)  c<^)(aj,  <^(a5)) ,  (22) 

of  the  medium’s  perturbation.  Here,  is  a  smooth  function  of  x,  some  (curved)  level  surfaces 
of  which  describe  a  family  of  interfaces.  The  gradient  of  the  perturbation  is  assumed  to  vary 
rapidly  normal  to  the  level  surfaces  of  <j>  and  smoothly  along  them,  implying  that 

Vajc^*)  =  (V®<^)  -h  terms  smoother  in  a; ,  (c^^^)'  =  .  (23) 

This  representation  extends  to  a  small  ball  around  any  point,  in  particular  the  image  point 
y,  say,  under  investigation.  The  derivative  is  understood  in  the  distributional  sense. 
Typically,  it  will  have  a  Dirac-distribution  type  behavior  across  any  geological  interface. 
Loosely,  the  derivative  can  be  interpreted  as  the  difference  in  medium  properties  across  a  level 
surface.  The  support  of  c^*l  is  denoted  by  V. 

In  the  background  medium  we  assume  that  the  ray  theory  of  the  previous  section  applies. 
In  the  fonvard  problem  the  background  c  and  the  perturbation  c^^^  are  both  known;  in  the 
inverse  problem  the  aim  is  to  reconstruct  given  c. 

We  wiU  omit  the  first  argument  of  in  the  remainder  of  this  paper,  and  take  only  the 
most  rapidly  varying  term  of  the  perturbation’s  derivatives  into  account  Imaging  reflectors 
or  interfaces  amounts  to  mapping  this  leading-order  behavior,  (c^^))'.  We  will  confirm  below 
that  our  inversion  procedure  does  this,  and  also  provides  estimates  of  angularly  dependent 
reflection  coefficients. 

4.  The  single  scattering  equation 

In  this  section,  we  introduce  the  Bom  approximation  representing  the  singly  scattered 
wavefield.  We  then  show  how  to  recast  this  volume  integral  representation  into  a  surface 
integral  for  the  response  to  the  most  singular  element  of  the  scattering  process,  namely, 
reflecting  the  discontinuities  of  the  medium  parameters. 

4.1.  Volume  integral  representation 

We  begin  the  analysis  with  the  volume-scattering  representation  of  the  ray-Bom  approximation 
for  the  scattered  displacement  field  for  the  {NM)  conversion.  Let  (r,  a)  e  dR  x  dS; 

ideally,  the  boundaries  dR,  dS  ~  =  unit  sphere)  are  closed  surfaces  surrounding  the 

heterogeneous  domain  V.  Then  (De  Hoop  et  al.  [14]) 

«W(r,.,i)  =  -jf  |(*)(,)(f)(.)x(»*)(*)  X  (24) 

(w(^^)(a.,a(^)(a.),a(")(a;)))^c(')(a;)(J"(t  -T(^^)(r,aj,s))  d*  , 


(25) 


where  N,M  €  {1,2,3}, 

=  p{x)  , 

contains  the  amplitudes, 

af)  =iv.'*'(f!*'7j'*>  +  jf'7j'^), 

describes  the  contrast-source  radiation  patterns,  and 


-I? 


represents  the  relative  medium  perturbation.  Here,  denotes  the  (local)  normal  speed  of 
mode  L  in  the  background  medium  averaged  over  all  phase  directions.  By  introducing  this 
convenient  scale,  we  have  made  the  quantities  a^i^\  and  dimensionless  since  the  7’s 
have  the  dimensions  of  [slowness]  and  cijki  has  the  dimension  of  [density]  x  [velocity]^.  The 
notation  « is  meant  to  emphasize  that  the  quantity  is  angle  independent,  which  is  important  for 
retaining  the  actual  medium  perturbation  from  and  the  GRT  inversion  to  be  applicable. 

In  the  micro-local  setting,  substituting  Eq.(22)  into  Eq.(24),  we  have 

t)  =  -//*)(r)ff  )(.)^''»*»)(*)  X 

(w(^J^)(aj,  d(^)(a5),  a(^)(*)))^  8"{t  -  r(a!))  d*  ,  (26) 

where,  for  convenience,  we  employ  the  shorthand  notation 

T{x)  =  s) ,  r(*)  =  X,  8)  ,  (27) 

see  Eq.(20).  To  make  use  of  the  properties  of  the  gradient  of  the  medium’s  perturbation 
(cf.  Eq.(23)),  we  will  partially  integrate  expression  (26).  Since 

(Va5r)(x)  5'\t  -  T{x))  =  -VxS\t  -  T{x)) ,  (28) 


we  have. 


8'\t  -  T{x))  = 


i>-VxS'it-T{x)), 


^  ^  "  (i>.v*r)(x)  ^ 

for  arbitrary  i>  G  as  long  as  i>  •  V*r  ^  0.  Hence, 

(w(^*)(*, d^^^(x), d^^^(*)))^  (c^‘^)'((/>(*))  X 
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Here,  the  approximation  arises  from  neglecting  lower  order  terms  as  in  Eq.(23),  as  well  as 
neglecting  derivatives  of  the  remaining  amplitude  in  Eq.(26),  compared  to  These  will 

produce  asymptotically  lower  order  contributions  to  the  wavefield. 

It  is  assumed  that  u  =  i>(x)  is  slowly  varying  in  space,  and  may  be  chosen  equal  to  the 
local  geological  dip. 


_  Va;</> 


(31) 


On  the  other  hand,  we  can  choose  v  =  v,  which  we  always  know.  In  the  inverse  scattering 
problem,  the  geological  dip  is  unknown  and  has  to  be  determined.  Below,  we  show  that  at 
stationarity  the  geological  and  migration  dips  must  be  parallel. 


4.2.  Surface  integral  representation 


Now,  we  show  how  to  recast  the  volume  integral  representation  (30)  into  an  integral  over 
surface  integrals  over  the  level  surfaces  of  <^.  To  this  end,  we  choose  curvi-linear  coordinates 
<r  =  o-{x),  <r  =  (<7^,  <73),  p  =  1,2,  such  that  the  are  coordinates  in  the  level  surfaces  of  (j) 
and  <73  is  the  local  coordinate  in  the  i/^-direction.  If  <73  represents  the  actual  value  of  the  level 
of  4>,  we  set  <73  =  L.  The  volume  form  is  given  by 

d*  =  ^  AL  dS(aj)  ,  dS(a;)  =  \d„,x  A  d„,xl  d<7id<72  ; 

the  nansformation  from  Cartesian  to  curvi-linear  coordinates  yields  the  Jacobian 


d{x) 

d{a) 


\da3X  •  {da,X  A  da^x)\ 


\d^,x  A  d<,^x\ 


(32) 


Now  write 


(cW)'(,^(*))  =  jf^(cC))'(£)i(*(*)  -  L)  di  , 


(33) 


Substituting  Eq.(33)  into  Eq.(30),  interchanging  the  order  of  integration,  and  using  the  property 
(cf.Eq.(32)) 

X...W«)-I)dx  =  /^^...^dE(x), 


0  -  -  /  f  /  X 

•/R 


(i>  •  r)  |V*^| 

=  -/d£[/  <<*)(r)|W(.)A'™>(»)  X 

•/R  L*'  4>^Lt 


we  obtain 
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{*>  •  t'4>) 

(i>.r) 


<x 


S\t-T{x))di:{x) 


(34) 


As  a  consequence  of  Eq.(34),  we  also  have 


a,«W(r, ,,t)=-[  ill f  x 

*/R  \.J  ih^lj 


(wW(x,dW(x),d<*)(aj)))^(c(*))'(L)  X  (35) 


(i>.r) 


The  singular  support  of  the  medium  perturbation,  <f>{x)  =  L,  and  the  isochrone  surfaces, 
r(x)  =  t,  are  illustrated  in  figure  3. 


geological  dip 


Figure  3.  Micro-local  medium  perturbation. 


5.  Reflection  and  transmission  coefficients 

To  identify  reflection  and  transmission  coefficients  in  Eq.(35),  we  first  extract  phase  velocities 
at  the  scattering  point  from  the  amplitudes, 

A„(x)  =  [v(^)(x)(V(^)(x))3]^^^  ; 


(36) 
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then,  Eq.(35)  can  be  written  in  the  form 

«(*!(*))  (i>  ■  ^,)(P  •  r)  u  S"(t 


in  which. 


(37) 


T{x)) 


dS(«)l 

|V*0|. 


dL, 


[K(iv)(y(M))3] '  |r|2  (i> .  uy 

represents  the  scattering  coefficient  for  the  ( N,  M)  conversion  at  x  with  (f){x)  =  L,  i.e., 
really  depends  on  <7^(05),  /x  =  1, 2.  Observe  that  the  scattering  coefficient  contains  a  function 
that  may  be  singular  on  the  level  surfaces  of  <i>.  In  fact,  the  integration  over  L  picks  up  the 
singular  support  of 


5.1.  Specular  reflection  and  transmission 

Now,  let  cri)  contain  a  step  function  in  L  (across  a  curved  interface  at  L  =  Lo).  Then 
(cri))'(Z,)  =  (Ac^^^)  5{L  —  Lo).  At  the  specular  point  for  given  i/^,  the  pair 
satisfies  Snell’s  law, 

^  •  (I  ^  ■  (I  -  ■  <39) 

[In  an  isotropic  medium  with  M  =  N  the  solution  is  simply  given  by 

representing  ordinary  reflection.]  At  the  specular  point,  the  migration  dip  coincides  with  the 
geological  dip,  1/  =  i/^. 

Substituting  the  solution,  dt[^\  of  Eq.(39)  into  the  scattering  matrix  Eq.(38),  yields  the 
linearised  reflection/transmission  coefficients 

&<'*>(.))  ■*(£  -  io)  =  «<'*'(■),  c«i*’(.))  •  (40) 

Here,  the  left  side  exploits  the  evaluation  of  at  specular  and  explicitly  notes  the 
distributional  character  of  the  right  side.  In  fact,  with  u  =  i/^,  Eq.(37)  is  equivalent  to 
the  Kirchhoff-Bom  approximation,  which  satisfies  the  principle  of  reciprocity.  For  notational 
convenience,  we  introduce 

Rf  *>(.,  cS%))  =  «('*>(.))  S(L  -  U)  (41) 

to  absorb  the  Dirac  distribution  in 
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5.2.  The  Kirchhoff  approximation 

Here,  we  show  how  to  go  from  the  linearised  Bom  representation  Eq.(35)  to  the  non-linear 
Kirchhoff  approximation.  For  general  reference,  we  introduce  the  relevant  scattering  and 


specular  angles:  in  addition  to  i/,  we  set 

cos  6  =  tj)  =  third  Euler  angle  .  (42) 

Then,  at  any  point  in  P,  we  have  a  mapping 

(43) 

If  1/  =  i/^,  the  ray  geometry  is  specular;  let  the  associated  specular  scattering  angles  be  defined 
as 

cos  6  =  cos  0  =  6^  =  6  —  0  .  (44) 

In  the  non-reciprocal,  ‘non-linear’  Kirchhoff  approximation,  the  scattering  matrix  is  simply 
replaced  by  iSa&jull  reflection/transmission  coefficients  at  specular,  cf.  Eq.(37)  with  u  =  v^, 

a,uW(r,.,i)  =  -£,[/  x  (45) 

*/J\  \,J  <p^Ij 

Rf  *>(*,&('»)(*))  (Ui  ■  r)  I,  3"((  -  T(x))  di  . 


For  wide-angle  (0)  scattering,  this  representation  is  certainly  more  adequate  than  Eq.(37). 
The  time  derivative  is  taken  to  pave  the  way  for  the  Radon  transform  inversion.  (In  three 
dimensions,  one  needs  the  second  derivative  of  the  Dirac  distribution.)  In  our  further  analysis, 
we  actually  employ  the  reciprocal  representation  Eq.(37)  in  which  is  replaced  by  the  full 

reflection/transmission  coefficient  at  specular.  Note  that  due  to  the  singular  function  contained 
in  cf.  Eq.(40),  the  integration  over  L  in  Eq.(45)  reduces  to  a  sum  over  scattering 

surfaces  or  interfaces. 

6.  Stationary  phase  analysis  of  the  direct  scattering  problem 

By  applying  stationary  phase  arguments,  the  integral  in  Eq.(37)  can  be  evaluated.  The  analysis 
confirms  the  consistency  with  asymptotic  ray  theory  in  configurations  with  a  family  of  surface 
scatterers. 

We  choose  rotated  Cartesian  coordinates  {x^,  z),  p  =  1, 2  in  the  neighborhood  of  a  yet- 
to-be  determined  specular  point  y{L)  e  {^  =  L),  such  that 

^  II  *'^7  -L  (46) 

and 
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see  figure  3.  The  function  T{y{L))  has  a  derivative  given  by 


dr 

dL 


|v*r| 


|V*^| 


»(i) 


(47) 


while,  by  the  implicit  function  theorem,  the  function  Ly{t)  satisfying 

T  (»(%{<)))  =  ‘ 


exists. 


Taylor  expansions  of  the  level  and  isochrone  surfaces  including  the  curvature  terms  yield 

0  =  “  l^sc^(y)l  ^  2^n^nu{,y)^v  i 


(48) 


T(aj)  =  T{y)d■\^xT{y)\zd■\x^J'^^{y)x^ 
(Summations  are  carried  out  over  fi,  u.)  Here, 


T„u{y)  = 


d^T 


y 


y 


(49) 


The  first  equality  in  (48)  amounts  to  the  representation  of  the  level  surface  {<f)  =  L}, 


Z 


2|V*<^|  ’ 

which  upon  substitution  in  the  second  equality  yields 


T{x)  =  T{y)  +  i|V*r(y)|a:^T^,(y)x, 


(50) 


and  z,  2f  in  the  same  level  surface.  Note  that  the  matrix  T  may  vary  with  the  level  L,  and 
can  be  negative  or  positive  definite,  or  indefinite.  The  case  of  vanishing  T,  leading  to  a 
caustic  analysis,  will  be  postponed  to  a  future  paper.  Otherwise,  for  t  near  T{y),  we  have  the 
intermediate  result 


-  r(x))  dE(x)  =  a?  -  t(*))  ds(x) 

(l  -  T{y)  -  5|Va.r(y)|  dn  di, 

_  27r  6*{t  —  T{y)) 

”  |r(»)lv'|det(T(»))|  ’ 
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using  tangent  plane  coordinates,  and  where 


(5  if  T  positive 

5*  ={  ns  if  T  indefinite 

—S  if  T  negative . 


(52) 


In  the  above  n  denotes  the  Hilbert  transform;  the  notation  *  indicates  an  action  on  the  time 
dependence. 

We  will  use  Eq.(51)  to  evaluate  the  integral  in  Eq.(34)  eventually.  First,  note  that  Eq.(51) 
implies 


27r 


|riv/|det(T)|  ("•r) 


l»(i) 


d-  [( -  r(»(£))] 


27r 


(i/  •  i/^) 


iriyjdit^  (i>-r) 

where,  for  any  function  f , 


viLyit)) 


(-V 

\dL)  ^  ’ 


(53) 


Ly{t) 


(S) 


-1 


=  /  f(£)d’l(-T(»{£))ldZ,. 


Substituting  Eq.(47)  into  Eq.(53),  and  using  the  result  in  Eq.(34)  implies 


\/|det(T)||r| 


(54) 


y{Ly{t))  ^  ^  \Ly{t) 


since  at  the  specular  point  we  have  u  =  and  u^  V  =  jr].  This  formula  is  an  extension  of 

the  convolutional  model  approximation  in  one-dimensional  space  to  three  dimensions. 

In  terms  of  the  reflection/transmission  coefficients,  we  have 


-  - 


27r 


yidSmiiri 


y  =  y{L) 

L  =  Ly{t) 


[note  that  the  *  relates  to  the  KMAH  index,  see  e.g.  Hdrmander  [22]].  To  verify  this  result, 
use  Eqs.(36),  (38),  (40),  and  (41)  in  (54). 
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7.  Inversion  based  on  the  GRT 


In  preparation  of  the  inverse  transformation,  we  introduce  the  scalar  quantity  on  the  diffraction 
surface, 

|W(r)  a(uW(r,  8,  y,  a)) 


dtu^^^\r,8,y)  = 

Then,  using  Eq.(37),  we  get 

dtu^^^\r,8,y) 


[v(J«)(„)(V'(A)(»))3] 


1/2 


-/[/  4“>(».6''*'(*).A<*V))  X 

«/R  L*f<f>=L 


(55) 


(i>  •  .  r)  I*  nny)  -  T{x)) 


dS(aj) 


dL.  (56) 


Here,  we  have  set  Au{x)/Au{y)  =  1,  which  is  its  value  at  the  dominant  critical  point  of  the 
integrand.  Furthermore,  we  will  exploit  the  expansion 

T{y)  -  T{x)  =  T{y)  •  (y  -*)  +  ...  ~  |r(y)|  (y  -*)•*/ ,  (57) 

which  describes  the  tangent  plane  to  the  isochron  at  y.  In  our  further  analysis,  we  will  make 
use  of  the  identity 

i"(|r(»)|  (y  -*)•.-)  =  |r(»)|-3y'((»  -*)■./). 

7.1.  Linearised  inversion 

The  basis  of  the  GRT  inversion  is  Gel’fand’s  plane  wave  expansion,  which  we  write  in  the 
form 

—  8Tr'^S{<f>(y)  —  L)  =  f  f  S{(f>(x)  —  L)S"{(y  —  x)  •  t/)  daj  du 
Js^Jv 


=  f  [  <^"((y  -x)-i^)  di/ 
Js^J<t>=L  |Vx0| 


(58) 


We  will  use  this  expansion  to  invert  Eq.(35).  The  inversion  is  accomplished  by  setting  up  a 
system  of  22  equations,  using  the  contraction  according  to  Eq.(55),  and  employing  a  and  d 
as  variables  of  integration,  i.e.. 


w(^*)(y,  a^^{y),  a^^\y))  s,  y) 


irr 


(l/  •  V4,) 


5(d,d) 


y 


dsdr 


y 


=  -/  dilf 

JH  U<I>=L 

w(^^)(y,  d<'^)(y),  d(^)(y))w(^^)(«,  d('^)(x),  6c^^\x)f{c^^)yiL)  x 


(i/-r) 

("  •  ^4') 

d{a,a) 

("  •  *'<!>) 

y 

(f-T) 

y 


^"{{y  —  x)'u)  dS(*)j di/d0dV> . 
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(59) 


E>efine  the  matrix 


di6t,a) 

d(v,e,rp) 


(60) 

dOd^ 


at  any  image  point  y.  Then,  employing  Eq.(58)  in  Eq.(59)  and  setting  i/  =  i/  (recall  that  u  has 
been  arbitrary  until  now),  yields 


(c">)'(^»))  |v*^|j,  =  i^y)  -  L)  \Vx<l>iy  H 

OTT^JdSxdR  ^ 


irr 


(l/  •  I/^) 


y 


d{a,6t) 

d{s,r) 


dsdr  . 

y 


X 


(61) 


Note  that  in  this  inversion  formula  the  actual  shape  of  the  level  surfaces  of  <}>  is  not  used; 
just  the  local  normal  at  the  image  point  plays  a  role.  However,  in  Section  9  it  will  be  argued 
that  the  inner  product,  (i/  *  i/^),  can  be  removed  from  the  formula  without  changing  the 
resolution  analysis  in  the  high-frequency  approximation.  The  inverse  of  A  is  understood  in 
the  generalised  sense.  We  will  denote  the  reconstruction  as  ((c^^^)'(^(y))  |Vaj«^|y). 

Through  Eq.(59),  we  have  the  freedom  to  carry  out  the  inversion  Eq.(61)  in  two  steps.  For 
a  given  image  point  y,  we  can  imagine  the  data  [(«,  r)  pairs]  to  be  sorted  into  common  {9,  ip) 
gathers.  The  variable  in  such  gathers  is  the  migration  dip  u;  formally,  we  denote  these  gathers 
by  {dS  X  dRY{y;  0,  ip).  The  integration  over  dip  is  then  carried  out  prior  to  the  integration 
over  scattering  angle  and  azimuth.  [In  practice,  shooting  the  rays  from  y  would  be  controlled 
by  a;  a  then  follows  from  (0, 0).  The  data  would  be  simply  taken  at  those  locations  where 
the  rays  intersect  dS  and  dR,  respectively.] 

To  control  the  illumination  of  the  image  point  at  a  given  dip,  we  introduce  the  partition  of 
unity,  {xj},  and 


((<:“>  W(y)) 

j  oTT  J dSXoR 


|r|^ 

(*'  •  ^<t>) 


Xj{8,r) 

y 


5(a,a) 

'dM 


dsdr  . 

y 


(62) 


Each  term  in  the  summation  represents  a  partial  reconstruction. 
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7.2.  Resolution  analysis 

Backsubstituting  the  Kirchho£f-Bom  scattering  formula  Eq.(3S)  with  1/  =  i/,  into  the  inversion 
formula  Eq.(62)  and  substituting  the  one-sided  Fourier  representation  of  the  Dirac  distribution, 

S"{T{y)  -  T[x))  =  -Re  -  (  exp[iu;  (T(y)  -  T{x))]uPAuj  , 

TT  JR+ 

yields  the  resolution  operator, 

((c<'>)'Wy)) 

JnJ<i>=L  |V*<p|® 

with  matrix  kernel  Tl, 

'7^(y ,  *)  =  53  exp[i$(y,  *,  0)]  a(y,  *,  0)  xj(0)d0  • 

Here,  0  =  (oj,  s,  r)  and 


(63) 


(64) 


i 


•  d0 


JdSxdR  In 


(65) 


•  |r(y)|^a;^da;  dadr  . 

fn+xdSxdR  JdSxdR  Jn^ 

The  resolution  operator  expresses  how  well  the  reconstruction  can  be  accomplished  within  the 
framework  of  the  linear  theory. 

The  phase  function  $  of  the  Fourier  integral  operator  with  kernel  Eq,(64)  is  simply  given 

by  . 

$(y,*,0)=a;(r(y)-r(*))  , 
while  the  amplitude  function  arises  as  the  matrix 

a(»,*,e)  H  ^§'»)(r,»)|f)(r,  *)  If  x 


(66) 


[(Ay(i/(r,  y,  a)))  ^  a^^^y),  Q:^^^(y)) 

(w(^^)(*,a(^)(a5),Q'(^)(aj)))^l  x 


■V(^)(y)(t>(^)(y))3‘ 

ir(y)|  (i/ •  i/^)aj  d{dt,6i) 

V(^)(*)(y(^)(*))3 

|r(*)|  (i/ •  i/^)y  a(a,r) 

y 


(67) 


In  anticipation  of  introducing  the  scattering  coefficients  (cf.  Eq.(38)),  we  introduce  the  one¬ 
dimensional  array  of  functions  a^,  satisfying 

an(y,  0)  a^^^(*)) 

=  a(»,  *,  0)  (c^^^)'(<;i>(*))  I Va5<?i>|aj  •  (68) 

We  will  interpret  0  in  the  spatial  Fourier  domain.  To  this  end,  we  carry  out  two  coordinate 
transformations.  First,  we  employ  the  ray-induced  mapping 

s  =  8{N,M,u,6,rp),  r  =  r{N,M,v,9,-^)  for  fixed  y  , 


(69) 
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with 


L 


d{dt,a) 


-L 


5(a,a) 


lasxdR  d{s,r) 

Thus,  at  y,  0  is  mapped  on  (w,  i/,  9,  V>),  and  we  set 

5(q!,q:) 


da  dr  ,  . .  ^  , 

y  Js»x!P  a(i/,0,V’) 


a(y,x,u,  I/, 0,11^) 


y 


di/  d^  d^  . 


=  a{y,x,u,8,r) 
y  d{i/,e,xp) 


y 


d{a,r) 

We  identify 

T^^^\r,x,s)  with  T^^^\x,t/,9,tp) . 

Second,  the  frequency  u  is  transformed  to  the  wavenumber  kr  according  to 

k,  =  u\T{y)\-,  (70) 

then 

•  •  •  ir(v)|V  do.  =  ■  t?  dtp .  (71) 

We  now  identify  the  wave  vector 

0'  =  krV  6  IR^  with  d0'  =  kf.  dfcr  di/  ,  (72) 

and  we  will  consider  {$,  t}))  as  parameters,  i.e.,  0  (w,  i/,  6, 0)  (0',  0, 0).  The  Jacobian 

of  the  latter  transformation  is  written  as  (cf.  Eq.(71)) 


d(0O 


5(aj,  u) 


=  h{y^ >  h{y^ *')  =  |r(y)|^ ; 


y 


formally,  also 


h  = 


det(r  duiT  5i/2r)|.  (73) 

The  inverse  transformation  to  frequency,  the  so-called  Stolt  mapping,  is  given  by 

„(e')  =  (74) 

since  also 

0'  =  u;r(y). 

We  identify 

$(y,aj,0)  with  $(y,®,0',^,0) . 

In  the  phase  space  with  coordinates  {x,  0'),  Eq.(64)  gives  rise  to  the  resolution  equation 

{(c<‘))Wy))|Vx^l,)  =  ReE/^^/X 

/  .  ^exp[i$(y,»,0',^,0)]a(y,a5,0',0,0)Xj(0',^,0)d0'  x 
(c('>)Wx))  |Vj,,M a,  di  dd  .  (75) 

|Va;0|x 
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For  X  near  y,  $  =  0'  •  (y  —  *)  H — .  The  integration  over  0'  represents  the  spatial  resolution, 
whereas  the  integration  over  6,  if;  primarily  represents  the  parameter  resolution  per  migration 
dip.  However,  note  that  the  parameter  resolution  couples  to  the  spatial  resolution,  since  0'  in 
general  depends  on  {0,  ip). 

Below,  we  will  constrain  xj  to  be  a  function  of  kr  =  |0'|,  ^  alone. 


7.3.  Stationary  phase  analysis  of  the  resolution  operator 


Inside  the  integral  over  w  in  Eqs.(64)-(65),  we  consider  u>  to  be  large.  Then,  we  apply  a 
four-dimensional  stationary  phase  analysis  with  respect  to  the  integrations  over  (oi,  02,  u)  G 
S{L)  X  5^,  cf.  Eq.(75).  The  range  S{L)  indicates  that  the  surfaces  are  contained  in  D. 

We  choose  polar  coordinates  on  the  i/-sphere. 


u  =  (sin  O'"  cos  Ip'",  sin  O''  sin  ip'',  cos  0'')  . 

We  extract  L  =  03  and  A:r  =  |0'|  from  the  set  of  phase  space  coordinates, 
(*,  0')  -)•  (<Ti,  <T2,  L,  kr,  9'', Ip’') ,  77  =  (cTj,  02, 0'',  Ip'') , 
resubstitute  kr  =  ci;|r(y)|,  and  set 

$(y,  X,  &,  e,ip)  =  u}  $'(y,  L,  t],  kr,  9,  ip) . 


Writing  the  coordinates  explicitly,  the  resolution  equation  (75)  takes  the  form 

((C'«)'W»))  |V.^|„)  =  R*  E  i  irJM  U. 
exp[iu^'{y,L,T),kr,9,ip)]a{y,L,r},kr,9,^)Xj{kr,9,ip)h{y,i/)  x  (76) 


where 


drf  =  dS(a5)  sin  9"' dO'' dip'' .  (77) 

For  given  (9,  ip),  the  phase  is  stationary  with  respect  to  the  variables  of  integration  77  if 

5^^$'  =  0  and  =  0  ;  (78) 

here, 

=  ,  (79) 

while 

(80) 
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The  stationary  points  are  denoted  as  and  induce  the  mapping  *(<7°,  L)  for  any  L. 

The  solution  of  the  first  equation  in  (78)  with  (79)  implies  that  the  stationary  migration  dip, 
must  be  parallel  to  the  geological  dip,  v^,  i.e., 

=  ±t'4> ;  (81) 

one  solution,  <t°,  of  the  second  equation  is  easy  to  identify: 

=  y ; 

at  this  value  one  expects  the  peak  contribution  to  the  resolved  medium  perturbation.  At 
aj(<T®,  L),  given  i/°,  {B,xj))  will  determine  the  stationary  values  of  (s,r).  We  denote  the 
stationary  points  by  rf  =  (<r°,  i/°(aj(<T°,  L))),  and  the  stationary  point  set  by  H®,  which 
contains  at  least  two  elements  (cf.  Eq.(81)). 

Applying  the  four-dimensional  stationary  phase  approximation  to  Eq.(76)  amounts  to 


exp  [i.  nv.  L,  K,  0,  ^)  +  sig(V,V,»')»] 


Xj{kr,e,xj})h{y,v^)  x 

(cW)  m 


(82) 

a;^du;dL  d^d^ 


in  the  absence  of  singularities.  [Singularities  require  a  separate  analysis,  which  we  will  discuss 
in  a  separate  paper.]  Here,  sig  denotes  the  number  of  positive  eigenvalues  minus  the  number 
of  negative  eigenvalues  of  a  matrix.  In  Eq.(82),  fcr  is  the  stretch  of  frequency  with  |r(  j/)  |°,  the 
norm  of  the  gradient  of  travel  time  in  case  the  migration  dip  is  stationary. 

In  terms  of  the  scattering  coefficients  (cf.  Eq.(68)),  for  *(<7°,  L)  near  y  ^  {(f>  =  £,}, 
expression  (82)  reduces  to 

((<:“’)'(«»))  |V*.^|y)  E  i  T,  L  L 

i/»  =  ±1/,  ■' 

Re  *'(»•  1°’ ''’)  +  •?  xAh, e, 

a/i(y,  L,  I)”,  tr,  '!>)  K,  P.  h(y,  i/°) 


|r(»)|»y|det(V,V,$')“l 

\d„^xNd„^x\ 


\^x<f>\ 


dL  dOdrl^ .  (83) 


Recognizing  the  bandlimited  Dirac  distribution, 

•  (y  -  *))  =  i  X 

TT 


(84) 
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/  exp  ftOJ  $'{y,  L,  7/°,  fcr,  0)  +  sig(  Xj(^r,  0)  , 

JR+  L  4  J 

we  find  that  the  reconstruction  amounts  to  a  weighted  integration  of  scattering  coefficients, 

((c''>)W»)) |v*^|j,) ^  T.  h'Ei  ,  L 

R* /i  jK  ■  (W  -  *))  X 

a«(y.  4 1?°.  <Jr,  0,  V*)  fijf ***(’?“.  <^r.  <i>)  My,  1^)  ^ 
|r(y)lVl*‘(^.^n»')“l 

|5»,x  dL  d0d0, 

where  we  have  accounted  for  the  fact  that  the  range  of  integration  over  9, 0  will  be  limited  by 
the  acquisition  geometry.  The  ranges  are  given  by  Ee  =  £0(1/°),  =  E^{i/°,  9). 

In  Section9,  we  will  evaluate  the  Hessian  det(V,,V,$')°andshowthatsig(V,V,$')°  =  0 

at  i/°  =  ±1/^.  Then,  in  case  the  partition  is  complete,  we  have 

E  ^  •(»-*))  =  jvbi 

~<S(0(y)-0(*)).  (86) 

Due  to  the  point  symmetry  of  the  slowness  surface  at  the  image  point  y,  we  can  replace  the 
summation  J2^o  _  \  by  the  substitution  i/°  =  u^. 

8.  GRT  inversion  of  Kirdihoff  data 


Now,  we  will  analyse  the  resolution  operator  from  a  different  perspective.  To  accommodate 
for  the  non-linear  reflection/transmission  coefficients,  we  substitute  our  Kitchhoff-like  approx¬ 
imation,  a  mixture  of  Eqs.(37)  and  (45),  into  the  inversion  formula  Eq.(61).  The  result  is  an 
equation  very  similar  to  Eq.(85).  We  obtain 

((c"’)'Wy))  IVx^ly)  E  5  (87) 

I/"  =  ±v^ 

where 

rl.“' Ra (s,-*))^ 

afl(y»  L,  rf,  fcr,  9, 0)  h,  9, 0)  /t(y,  u°)  ^ 

|r(y)lo^|det(V,V,$01 

la^.a;  .  (88) 


22 


Componentwise,  the  reflection/transmission  coefficients  can  be  identified,  viz.,  by  undoing  the 
multiplications  by  a^. 

9.  Imagiiig  reflectivity 

9.1,  The  Hessian 

In  this  section  we  will  evaluate  the  Hessian  of  at  a  stationary  point.  First,  note  that 
=  0  at  the  stationary  point  x  =  y  (then  $  =  0).  Hence, 

det  ( V, =  [det  ( V,/ Vtr$')“] "  •  (89) 

On  the  other  hand,  note  that 

•  A.x) .  (90) 

This  matrix  can  be  written  in  the  form 

(  I  tr  I  )  (  Y  Y  )  ’ 

hence,  also. 


-  -  ] 

(I  1  1  Y 

|r|  detPi/F)  •  (da-x)]  =  det 

-  d^r  - 

a.,®  s„* 

[-  T  -J 

1  1  1  1  /J 

(91) 

In  this  expression,  the  first  matrix  on  the  right-hand  side  can  be  identified  as  h  (see  Eq.(73)) 
and  the  second  one  with  the  Jacobian  of  the  coordinate  transformation  Eq.(32).  Hence,  using 
Eq.(89),  we  arrive  at  the  identity 

iri"  =  *(.,./“)  la.,*  A  .  (92) 

In  view  of  Eq.(89),  the  positive  and  negative  eigenvalues  must  come  in  pairs,  so  that 
sig  ( V,  V,$')°  must  be  equal  to  0  or  4.  On  the  other  hand,  we  have  intrinsically  assumed  that 
the  gradient  of  two-way  travel  does  not  vanish,  so  that  /i  ^  0.  In  accordance  with  Eq.(92), 
hence,  the  determinant  caimot  vanish.  This  implies  that  under  continuous  deformations  of  the 
interface  and  ray  geometries,  the  signature  of  the  Hessian  cannot  change  from  0  to  4  or  vice 
versa  (this  would  require  an  eigenvalue  to  become  zero).  In  the  case  of  flat  level  surfaces,  it 
can  be  shown  that  the  signature  equals  zero,  which  now  implies  that 

sig(V,V,$0°  =  O 


for  any 
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9.2.  The  modified  GRT  inversion 


Substituting  Eq.(92)  into  Eq.(88)  now  yields 

^(NM)  =  ^  Re  /a,j(i/°  •  {y  -  x)) 


x{crl,L) 

aji(y>  L,  ri°,  kr,  0, 0)  kr,  0,  ip)  , 


(93) 


which  is  the  proper  interpretation  of  the  set  of  images  created  by  the  GRT  inversion  formula 
(61).  To  arrive  at  this  expression,  we  could  use  that  at  the  stationary  dip  i/°  =  1/4,,  and  we 
could  set  {1/  •  1/4,)  =  1  in  Eq.(61)  to  begin  with.  Thus,  a  priori  knowledge  about  the  geological 
dip  is  not  required. 

In  the  stationary  phase  approximation,  we  can  rewrite  inversion  formula  (61)  with  Eq.(87) 


as 


Jji 


Here,  we  employ  the  mappings  dehned  in  Eq.(69), 

To  extract  the  reflection  coefficient  from  one  has  to  estimate  the  a/?,  which  are 

functions  of  the  stationary  dip,  scattering  angle  and  azimuth.  In  the  inversion  procedure,  we 
control  the  scattering  angle  and  azimuth;  in  principle,  we  can  estimate  the  geological  dip  from 
any  of  the  images.  )^th  this  estimate,  the  can  be  evaluated.  Now,  note  that  (93)  comprises 
a  system  of  equations;  from  each  equation,  in  principle,  the  reflection  coefficient  at  specular 
can  be  determined.  This  redundancy  can  be  employed  to  verify  or  improve  the  estimate  of  the 
stationary  dip.  The  stationary  dip  also  appears  in  the  spectrum  of  the  medium’s  perturbation; 
this  is  discussed  in  Appendix  A. 

On  the  other  hand,  by  virtue  of  the  stationary  phase  approximation,  we  can  remove  the 
AVA  inversion  nested  in  formula  (61):  consider  the  procedure 


1  r  \T{y)\^  y) 

Sn^Js^  [VW(y)(V(M)(y))3]'/" 


di/ 


5 


then  a/i  must  be  replaced  by  the  scalar  quantity 


SiR{y,x,e)  ->■ 


^u(y) 


Its  diagonal  is  given  by 


(95) 


(96) 


^R{y,y,Q)^  |r(y)r  ^ 


(97) 
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By  producing  images  both  with  m  =  4  and  m  =  3,  |r(y)|  at  stationary  can  be  imaged  as  well, 
viz.,  from  the  ratio  of  the  images.  Then  the  stationary  dip  is  not  required  to  find  the  reflection 
coefficient 

10.  Discussion 

We  have  shown,  by  carrying  out  a  stationary  phase  resolution  analysis,  that  it  is  feasible  to  ex¬ 
tract  information  about  the  angular  dependent  reflection/transmission  coefficients  from  a  GRT- 
based  migration/inversion.  We  did  not  have  to  linearise  nor  expand  the  coefficients.  In  fact 
the  outcome  of  the  resolution  analysis  is  a  multiple  set  of  images  for  the  reflection/transmission 
coefficients  for  the  available  range  of  specular  scattering  angles.  Any  type  of  AVA  analysis  can 
then  be  applied  to  interpret  those  images.  In  the  derivation  we  have  made  use  of  the  fact  that 
the  surface  integral  representations  are  linear  in  the  scattering  coefficients;  these  coefficients 
reduce,  at  specular,  to  the  reflection/transmission  coefficients. 

The  GRT  approach  employs  a  somewhat  unusual  input  of  data,  viz.,  via  common  (d,  0) 
gathers.  The  inversion  formula  reduces  to  a  two-dimensional  integration  over  migration  dip. 
The  (0,  V’)  sorting,  however,  varies  with  the  image  point.  It  bears  resemblance  with  the  sorting 
in  common  offset,  though.  The  use  of  such  a  sorting,  however,  necessitates  the  calculation  of 
an  additional  Jacobian. 
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Appendix  A.  The  spectrum  of  the  medium  perturbation 

In  principle,  the  geological  dip  can  be  directly  estimated  from  an  image.  However,  the 
geological  dip  is  also  hidden  in  the  spectrum  of  die  medium  perturbation  Eq.(33). 

Setting  (cf.  Eq.(46)) 

^  II  ^4-y  {®/*}  -L 

at  y  G  D,  the  medium  perturbation  spectrum  is  of  the  form  (cf.  Eq.(33)) 

d^i)(fc)=  /  dI(c(^)y(L) 

•/R 

j^S{<l){x)  —  L)  exp  [-i  (k^Xf^  4-  k2{zy{L)  -f  z))]  d*  ,  (Al) 
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where  in  fact  fc  —  0'  =  k^u.  Near  the  level  surface  (^  =  L,  we  employ  the  expansion 

^(aj)  =  L  +  \Vx<i>\z  +  ;  (A2) 

we  implicitly  assume  that  the  integral  over  L  is  windowed  around  y  such  that  in  the  window 
or  neighboriiood  may  depend  on  L  but  does  not.  At  the  level  surface  on  which  y  lies, 
we  have  zy  =  0.  Then, 


5^cW(fe)  ^  j^dI(c^^)y(I)exp[-iA;2Zy(L)] 

jm  H"  ^  “h  dxi  dx2  dz 

=  ^  exp[-iA:^^y(L)] 

/  exD  \-i(k  X  -k 
_  f  27r  exp  [7rtsig(A:^<^^^)/4] 

Jr 


.■>1 


ic^^^y{L)exp[-ik^zy{L)]  dL  ,  (A3) 

where  sig,  as  in  the  main  text,  represents  the  sum  of  signs  (±1)  of  eigenvalues  of  Note 

that  kii  =  0  corresponds  with  the  geological  dip  direction. 

If  the  level  surfaces  of  <f>  were  flat  and  t/^  =  fixed,  we  would  get 

9^c(^)(fc)  =  8{k  -  (k  •  (c(^))'(jk  •  i/^) 

^  ■  *'^^*'*^  (c(i))'(A:r(*/  •  i/^))  , 

which,  in  the  Radon  domain,  implies 


J^d^c^^\4>{x))  S"{y  •  1/  -  a;  •  I/)  d*  (A4) 

=  --^(i/-(i/-*/^)i/^)Re^^(c(i))'(fcr(i/-*/^))  exp(ifcry’)  dfc, 

=  8{u  —  {u  •  u'-)u^)  R«^^(c(^))'(A;r)  exp(ifcr¥?)  dfc, 


(pz=y.u 


ip  =  y.u 
1  • 


This  formula  shows  that  the  GRT  algorithm  can  reveal  the  geological  dip  explicitly.  We 
assumed  a  proper  coordinate  system  on  5^,  such  that 


-  (i/  •  di/  =  1 . 
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Figure  captions 


Figure  1.  Source-receiver  ray  geometry  (5^  =  unit  sphere). 


Figure  2.  Source-receiver  ray  polarizations. 


Figure  3.  Micro-local  medium  perturbation. 


